Computational analysis of 5-fluorouracil anti-tumor activity in colon cancer using a mechanistic pharmacokinetic/pharmacodynamic model

5-Fluorouracil (5-FU) is a standard chemotherapeutic agent to treat solid cancers such as breast, colon, head, and neck. Computational modeling plays an essential role in predicting the outcome of chemotherapy and developing optimal dosing strategies. We developed an integrated mechanistic pharmacokinetics/pharmacodynamics (PK/PD) model examining the influence of 5-FU, as an S-phase specific double-strand break (DSB)-inducing agent, on tumor proliferation. The proposed mechanistic PK/PD model simulates the dynamics of critical intermediate components and provides the accurate tumor response prediction. The integrated model is composed of PK, cellular, and tumor growth inhibition (TGI) sub-models, quantitatively capturing the essential drug-related physiological processes. In the cellular model, thymidylate synthase (TS) inhibition, resultant deoxynucleoside triphosphate (dNTP) pool imbalance, and DSB induction are considered, as well as 5-FU incorporation into RNA and DNA. The amount of 5-FU anabolites and DSBs were modeled to drive the kinetics of the pharmacological tumor response. Model parameters were estimated by fitting to literature data. Our simulation results successfully describe the kinetics of the intermediates regulating the 5-FU cytotoxic events and the pattern of tumor suppression. The comprehensive model simulated the tumor volume change under various dose regimens, and its generalizability was attested by comparing it with literature data. The potential causes of the tumor resistance to 5-FU are also investigated through Monte Carlo analysis. The simulation of various dosage regimens helps quantify the relationship between treatment protocols and chemotherapy potency, which will lead to the development of efficacy optimization.


Introduction
In this study, we developed a novel mechanistic PK/PD model that integrates a detailed cellular model with the classical PK/PD models to analyze the treatment outcome of the chemotherapeutic agent 5-fluorouracil (5-FU) in colon cancer. The modeling framework enables the quantitative analysis of the drug metabolism together with the downstream antiproliferative events by providing process models for the distribution of 5-FU from its entry into the body to interstitial fluid outside the tumor cells and to its cellular uptake. In the cell model, the connection of cellular level molecular events that lead to tumor suppression are modeled in detail. Over the past decades, pharmacokinetics/pharmacodynamics (PK/PD) models have been widely used in predicting therapeutic outcomes through the evaluation of the causal concentration-effect relationships [1]. Specifically, the PK model concentrates on drug disposition in plasma and/or tissues and describes drug exposure stages. In contrast, the PD model examines the pharmaceutical responses (e.g., tumor volume, white blood cell counts) following the drug treatment. The structures of the PK models differ based on the type of the drug, route of administration (i.e., oral absorption, intravenous injection, intravenous infusion), as well as the availability of clinical data. Commonly used compartment PK models in the literature are one compartment, two-compartment and three-compartment models. The one-compartment PK model considers the drug distribution in the whole organism as homogeneous [2]. The two-compartment PK model represents the transfer of substances between central (richly perfused organs) and peripheral regions (poorly perfused organs) of the human body [3,4]. One central and two peripheral compartments constitute three-compartment PK models [5,6]. The multicompartment pharmacokinetics models are preferred over the one-compartment model, because they can accurately represent the plasma concentration-time profile [7]. The pharmacokinetic profile in the central or peripheral compartment is the input of the PD model. The drug transport across the cell membrane is modeled as the irreversible or reversible exchange with the tumor compartment [8]. The PD model for tumor response comprises the kinetics of the control group and the treated group. The population tumor growth kinetics in the absence of the drug intervention is described by ordinary differential equations, including a growth term and a death term. Exponential, Mendelsohn, linear, logistic, Gompertz, and Bertalanffy are the common options for simulating the normal tumor growth patterns [9]. Specifically, exponential, Mendelsohn, linear, and surface growth terms predict that tumors will continue growing without an upper bound. Logistic, Gompertz, and Bertalanffy models predict that tumors will grow to some maximum size and reach a stable equilibrium at that point, based on the finding that the tumor growth is subject to the limited nutrition supply and extracellular cues [9,10]. Under the effect of tumor therapy, the tumor growth is altered and deviated from its baseline level. As a key player in PD model, the exposure-effect relationship (E drug ) driving the drug-induced tumor kinetics is modeled as a function of drug concentration, which can be linear, log-linear and sigmoid [11,12]. The involvement of drug effect in tumor growth dynamics can be modeled in various ways, depending on the biological question researchers propose to address, the assumptions on which the model construction is based, and the availability of experimental data. In general, the pharmacodynamic models take one or more of the forms of the hypothetical effect compartment model, indirect response model and signal transduction model, irreversible effect model, and complex mechanism-based model [13]. The interpretations of these PD models are provided as follows. The hypothetical effect compartment reproduces the true site of drug action and is applied between the PK and PD model only when the slow onset of observed effects can be explained by the rate-limiting drug distribution to the target site [14]. The indirect response model captures the nonlinear exposure-effect relationship and reflects how the reversible drug action alters the turnover of the drug response variables. Specifically, the sigmoid E max,damage can be added to the production and loss term to represent either inhibitory or stimulatory effects of drug disposition on the production or the loss of the responses [15,16]. The PD model encompassing the signaling pathways incorporates the time-dependent intermediate steps connecting the drug-receptor complex and pharmacodynamic response. The time elapsed from drug binding to the onset of the measurable responses is modeled by a transit compartment approach [17,18]. The irreversible effect model is similar to the indirect response model in certain contexts where the tumor growth is incorporated. But this type of PD model is only applicable to the compounds with irreversible anti-tumor effects and does not alter the turnover rate of a biomarker. The mechanism-based pharmacodynamic models extend from the indirect response model and consider the multiple drug-related mechanistic events represented by the ordinary or partial differential equations defined for the model variables [19]. The underlying physiological processes also warrant the additional modeling components to capture the biological phenomenon.
In recent years, researchers have shown great interest in employing PK/PD models for the prediction of tumor responses to 5-FU. Several conventional PK/PD models in the context of 5-FU treatment have been proposed, in which the transit compartment model is employed to describe the progressive effect on damaged tumor cells [20][21][22]. Although researchers have gained desirable analysis outcomes through the PK/PD models, the physiological effects under tumor treatment are barely captured. To delve into the complex cellular responses towards chemotherapy and to realize precise prediction, the mechanistic details need to be incorporated in the model design phase. In the literature, semi-mechanistic PK/PD models have been built for 5-FU in recent years, which employ both biological complexity and predictive power of a typical empirical model. In the model developed by Arshad et al. [23], three transit compartments describing blood cell maturation were used to explore the relationship between 5-FU exposure and myelosuppression. In a physiologically-based PK model [24], the activities of the metabolic enzymes regulating the metabolism of Capecitabine to 5-FU are considered to pursue a better prediction of the therapy efficacy. The study [25] aiming for characterizing the PK/PD relationship for LY2835219 developed an integrated model incorporating CDK4/6 inhibition, cell-cycle arrest, and tumor growth inhibition. Altinok et al. proposed a cellular automaton model coupled with a PK model to investigate the effects of circadian rhythms on the therapeutic outcome of 5-FU [26]. Recognizing the benefits brought about by this type of computational model, we developed a mechanistic PK/PD model which encompasses multilevel biological phenomenon to capture the drug actions of 5-FU in the tumor cells.
5-FU exhibits its cytotoxicity through a metabolic pathway that has been widely studied since its approval. 80%-90% of 5-FU is further decomposed into simple molecules through catabolism taking place in the liver by dihydro-pyrimidine dehydrogenase (DPD), while its active forms of metabolites are produced through anabolism in tumor cells [27]. It is 5-FU anabolites that intervene in the normal cellular functions and exhibit deleterious effects on tumor growth. Fluoro-deoxyuridinetriphosphate (FdUTP) and fluorouridine triphosphate (FUTP) are the products of sequential phosphorylations initially catalyzed by thymidine phosphorylase (dThdPase) and uridine phosphorylase (UrdPase) or orotate phosphoribosyl transferase (OPRT), respectively [28,29]. They are ultimately incorporated in RNA and DNA through the action of DNA polymerase, causing genomic instability. The conversion of 5-FU to fluorodeoxyuridine (FUDR) is followed by the production of fluorodeoxyuridylate (FdUMP) catalyzed by thymidine kinase (TK). FdUMP binds to thymidylate synthase (TS) at its nucleotide site in competition with deoxyuridine monophosphate (dUMP) and forms a ternary complex with a folate analog 5,10-methylenetetrahydrofolate (CH2THF) as a co-substrate, leading to reversible inactivation of TS enzymatic activity and inhibition of thymidine production [30,31].
Accumulating evidence indicates that TS inhibition is the determinant of 5-FU's antitumor effects and its activity has been demonstrated to be the predictor of tumor response to 5-FU treatment [31][32][33]. Hence, in our model, the dNTP pool imbalance caused by TS inhibition is attributed to the key factor in the formation of deleterious DNA fragmentation. The intactness of the components in the dNTP pool is required for DNA replication and DNA repair machinery to perform with high fidelity [34]. The thymidine deprivation and resulting dUMP accumulation can give rise to the perturbation of the biosynthesis of the other deoxynucleotides (dATP, dGTP, dCTP) via multiple feedback mechanisms. Interruption of biosynthesis of DNA precursors caused by the deficiency of the enzymes involved in nucleotide metabolism can give rise to futile DNA damage repair and accumulated DSBs [35]. (F)dUTP insertions into DNA facilitated by dUMP accumulation are recognized and excised by uracil DNA glycosylase (UNG), the first responder of the BER pathway [36]. But the uracil bases are incorporated repeatedly, causing the futile cycling repair by BER [36]. The accumulation of BER intermediates such as AP sites and single-strand breaks (SSBs) and persistent breakages triggers the formation of DSBs, resulting in cell cycle arrest, activation of homologous recombination repair(HRR), and apoptosis pathways. Indeed, some studies have established the relationship between the dNTP pool imbalance and 5-FU cytotoxicity [32,37].
In this study, we constructed a detailed mechanistic model linking the concrete drug actions of 5-FU to the prediction of pharmacodynamic response of colon cancer cells. The cellular model focuses on 5-FU anabolism, RNA and DNA misincorporation, TS inhibition leading to contamination of the dNTP pool, and induction of DSBs when the cell population proceeds through the S phase, all of which are believed to be essential biological events upon the administration of 5-FU. In addition, we used the γ-H2AX foci count to approximate DSBs, γ-H2AX is the canonical biomarker of DSBs, and its qualitative consistency with DSBs has been well established in the literature [38]. The tumor growth inhibition (TGI) model is built based on the combination of the life span model and the indirect response model, and the dynamics of the normal tumor growth is captured by the logistic growth model involving carrying capacity. Considering the suppressive effect of drug on tumor growth, 5-FU anabolites are modeled to inhibit the tumor growth and 5-FU-induced DSB formation is modeled accelerate the conversion of proliferating cells to damaged cells, which quantitatively captures the observations that the induction of DSBs by 5-FU and 5-FU anabolites are candidate predictive markers of 5-FU treatment efficacy [39,40]. In addition, the TGI model is constructed and validated using literature tumor volume data observed from the tumor xenografts administered with different dosage regimens. Another advantage of our mechanistic model is that the simulation-based resistance analysis can be conducted to provide the computational assessment of the potential causes of tumor resistance to 5-FU treatment. We not only examined the role of individual reactions in causing tumor resistance to 5-FU but also took into account the interactions among multiple factors. By comparing the time courses corresponding to these simulation scenarios and the time course of the tumor volume measured in 5-FU resistant cell line, we computationally identify which events contribute significantly to the drug resistance. Furthermore, we carried out global sensitivity analysis by applying the Morris method to screen out the subsets of the parameters with considerable effects on model outputs and then applying Sobol analysis that allows for the analysis of parameter interaction effects.

Materials and methods
Our mechanism-based model is composed of PK, cellular, and tumor growth inhibition (TGI) sub-models, the last two of which constitute the PD model. The model diagram is shown in Fig 1. The PK model is used to compute the concentration-time profile, and the TGI model focuses on the kinetics of tumor growth post-treatment. Our cellular sub-model captures physiological processes at the molecular level and connects the drug distribution in plasma and pharmacological response towards the drug action.

PK model
The PK model is described by a two-compartment plasma-peripheral model with Michaelis-Menten clearance [8,[41][42][43]. The central compartment models the drug distribution in the central plasma. All the remaining tissue space is captured by a hypothetical compartment, referred to as the peripheral compartment so that the exchange of drug substances between central blood (compartment 1) and auxiliary region of human body (compartment 2) given by a linear kinetics, and non-linear elimination collectively describe the drug fate in the circulatory system.
where Q 12 is volume of 5-FU cleared inter-compartmentally between central and peripheral space per minute. V 1 and V 2 are the volumes of the central and peripheral compartments, respectively. V max 1 and K m1 describe the process of elimination of 5-FU from central plasma.

Cellular model
The cellular model is built using a compartmental modeling approach, which is one of the common techniques utilized to model biochemical pathways. In the previous model involving similar molecular complexity developed by Wolf et al [44], the drug transfer from interstitial fluid to intra-cellular space and the emergence of products of drug action are considered. In our work, we construct a model including a more detailed description of physiological events associated with 5FU cytotoxicity. The cellular model considers the 5-FU concentration in interstitial fluid (compartment 3), formation of 5-FU anabolites (compartment 5), insertion of FUTP, FdUTP into RNA and DNA (compartment 6 and 7) binding of FdUMP to TS (compartment 9), changes in the components of the dNTP pool caused by reduced TS catalytic activity (compartment 10) and DSB induction as a response to imbalanced dNTP pool (compartment 11). For the purpose of modeling, the compounds are assumed to dissolve in a

PLOS COMPUTATIONAL BIOLOGY
Computational analysis of 5-fluorouracil anti-tumor activity spatially homogenous environment and be confined in the compartments, in which certain reactions with the respective compounds as reactants would take place. The rate of change of these reactions are represented by the influx and efflux of corresponding compartments. The dynamics of those species are subject to a group of differential equations described by the law of mass action and the Michaelis-Menten kinetics. The time evolution of physiological intermediates, as the output of the cellular model, is involved in the determination of drug effect measures, driving the ultimate pharmacodynamic responses. 5-FU in interstitial fluid and its anabolism. The differential equations describing the evolution of 5-FU in interstitial fluid (A 3 (t), pmol�mg −1 ), intra-cellular 5-FU (A 4 (t), pmol�mg −1 ) and 5-FU anabolites(A 5 (t), pmol�mg −1 ) are formulated in Eqs (4), (5) and (6). The kinetics of 5-FU in interstitial fluid (A 3 (t), pmol�mg −1 ) is governed by unidirectional transmission from central plasma [45,46], first-order elimination from the extracellular domain (k 03 ) and update and removal of 5-FU molecule across the tumor cell membrane [46]. The influx and efflux transporters anchored on the cell membrane contribute to uptake and clearance of 5-FU drug molecules [47]. Considering that the drug diffusion into the target cell is limited by the maximum number of available uptake transporters, the Michaelis kinetics, characterized by the saturation curve, is applied to describe the process of absorption of 5-FU by tumor cells. It is difficult to detect the individual low molecular-weight 5-FU anabolites (i.e. fluoronucleotides and fluoronucleosides, denoted as FNUC) accurately despite using high precision techniques [48]. Therefore, FNUCs are collectively represented by variable A 5 (t) and are modeled to be uniformly distributed and non-diffusing in one compartment. The rate of formation of A 5 (t) is through 5-FU anabolism. The decomposition of the TS-FdUMP complex and removal of 5-FU from RNA account for the replenishment of A 5 (t). The decrease of amount of A 5 (t) is regulated by the incorporation of FdUTP and FUTP into DNA and RNA and attachment of FdUMP on TS.
The time-dependent levels of incorporated FU residues into RNA(A 6 (t), FU-RNA, and DNA((A 7 (t), FU-DNA) are the solution of Eqs (7) and (8). The influxes into compartment 6 and compartment 7 quantitatively capture the incorporation of 5-FU into RNA and DNA, receptively, the rates of which are driven by the kinetics of 5-FU anabolites (A 5 (t)). The efflux terms in Eqs (7) and (8) describe the efficiency of excision of uracil bases by cell-intrinsic repair machineries. Different from the linear kinetics for FU-RNA, Michaelis-Menten kinetics is applied to FU-DNA since FdUTP insertion into DNA and restoring FU-containing DNA is regulated by protein enzymes. According to the experimental data regarding 5-FU insertion into genomes, the process of removing genomic 5-FU is not instantaneous but rather lasts for days. It is also shown in the literature that the initiation of cell cycle checkpoint and futile cycle DNA repair are the two main causes of delayed DNA repair process. The persistence of uracil analogs in genomes challenges genomic stability, initiating the cell cycle checkpoint and preventing cell cycle progression. The cellular repair machineries are precipitously assembled. However, the futile cycling repair, the causes of which have been discussed at length in the introduction, leads to the accumulation of unrepaired mismatched bases and impeded operation of DNA repair proteins. In conformity with the underlying mechanism, the rates of loss of FU-RNA and FU-DNA are the function of past states of FU-RNA and FU-DNA, realized by two delaying parameters (T d,RNA , T d,DNA ). Specifically, for t � T d,RNA (T d,DNA ), the loss of genomic 5-FU in RNA(DNA) is modeled as the kinetics without delay, while, for t > T d,RNA (T d, RNA(A 6 (t), pmol�mg −1 ) and FU-DNA (A 7 (t), pmol�mg −1 ) are shown below.
In which C p (t) denotes concentration of 5-FU in central plasma, which is solely defined by the PK model, Q 31 is the volume of 5-FU cleared from central plasma to extracellular space per unit time, V max,influx along with K m,influx describes the saturable uptake, V max,efflux and K m,efflux describe the saturable removal, k 03 is the first order rate of eliminating 5-FU from interstitial fluid, V max,54 is the maximum rate of 5-FU anabolism, and K m54 is the corresponding half saturation value. k 56 and k 65 are the first order rates of incorporating FUTP into and removing FUTP from RNA. k 06 and T d,RNA as well as k 07 and T d,DNA , describe the delayed excision of incorporated FU residues from RNA and DNA. γ lag,RNA and γ lag,DNA are used to adjust curve shape. The DNA polymerase with limited catalytic capability controls the insertion of FdUTP into DNA, expressed by parameters V max,75 and K m, 75 . V max,7 and K m,7 represent the rate of removal of FdUTP from DNA by UNG. Initial conditions for each compartment are zero. The units of incorporation of 5-FU into DNA and RNA are fmol� μg −1 DNA and pmol� μg −1 RNA, as shown in the data source literature, which are converted to pmol� mg −1 tissue through the concentration of DNA and RNA(0.2 μg/mg tissue and 2 μg/mg tissue respectively) [49].
The next dynamics we need to model is the binding of FdUMP to TS, which is considered to be the main contributor to 5-FU cytotoxicity. The study on analyzing the protein structure of E.coli TS [50] found that the methylation of dUMP to dTMP by TS with CH2THF as the methyl donor is a sequential process in which the dUMP binds at the active site of TS before CH2THF does. However, from the quantitative perspective, the differences in rate between the nucleotide/folate-ordered and random binding mechanism of reaction at TS is minimal [51]. Thus, for the sake of simplicity, we assume that the reaction involving dUMP is subject to rapid-equilibrium random mechanism instead of the ordered bi-substrate reaction. Since FdUMP and dUMP share the common biochemical properties in terms of covalent binding to TS [52], the pattern of interaction between FdUMP and TS is also in the manner of rapid equilibrium, demonstrated by the first term in Eq (9). It is also assumed that the free TS is saturated with CH2FH4.
To model the rate of TS-FdUMP complex formation, the quantity of TS enzyme over time is also considered, as shown in Eq (10). It is a consistent finding that the TS expression increases in response to 5-FU treatment [53][54][55]. TS expression at both protein and mRNA levels are regulated by multiple factors such as activation of cyclin/CDK complexes and its auto-regulatory feedback [56]. The increase in TS level after 5-FU exposure, with contributing factors still indefinite, may be related to gene amplification and disrupted synthesis regulation [56]. Therefore, it is more appropriate to estimate the total empirically content of TS by a mono-exponential equation.
Free FdUMP binding sites (TS f ) is expressed in terms of the all available FdUMP binding sites (TS total ), and FdUMP-bound TS enzyme (TS-FdUMP). In the literature that we use as the data source, the experimental measurements of TS total are carried out after and TS f are carried out before the dissociation of the ternary complex respectively [57]. Pretreatment TS f level is set as initial level of TS total . The amounts of TS total and TS f upon exposure to 5-FU are determined by Eqs (10) and (11).
The kinetics of TS-FdUMP complex(A 9 (t), pmol�mg −1 ) is expressed by Eq (9). Eqs (10) and (11) are the analytical expressions for TS f and TS total . dA 9 ðtÞ dt ¼ k 95 ðTS total À ½TS À dUMP�ÞA 5 ðtÞ À k 59 A 9 ðtÞ À k 09 A 9 ðtÞ ð9Þ In which the rate of enzyme-substrate complex formation is denoted as k 95 . k 59 represents the complex dissociation rate. The complex dissociation is accompanied by the release of FdUMP and the recovery of free TS. k 09 represents the first-order process of complex degradation. The production rate of the TS-FdUMP complex is proportional to (TS total − [TS-dUMP]), implicating that the competition between dUMP and FdUMP on the TS nucleotide binding site impedes the attachment of FdUMP to TS. α is the ratio of increased TS synthesis rate in the presence of drug to the normal synthesis rate of TS, k d is the ratio of degradation rate measured in treated cells and normal degradation rate constant, TS 0 is the level of TS protein prior to 5-FU administration.
It is assumed that there is an instantaneous equilibrium relationship between free dUMP and TS-dUMP-folate ternary complex. Therefore, the complex [TS-dUMP] is given by Eq (12).
where K dUMP is the dissociation constant for the enzyme-substrate complex and is the ratio of reverse and forward rate constants. It is more convenient to express complex [TS-dUMP] as a function of free TS and the total dUMP(dUMP total ) that is the sum of free dUMP and [TS-dUMP] complex, since the quantities of free TS and total dUMP are variables considered in the cellular model. By solving free dUMP out of total dUMP and [TS-dUMP], we can get the new expression for variable dUMP f that is dUMP total � . Following this, [TS-dUMP] can be further expressed in terms of dUMP total and TS f , as shown in Eq (13).
where G 0 is the zero-order rate constant of dUMP synthesis. k 08 is the first-order degradation rate of dUMP. r TS represents the rate of conversion of dUMP to dTMP catalyzed by TS, which is responsible for the loss of total dUMP. The expression of r TS is shown in Eq (15) and is linearly dependent on the amount of dUMP/TS/folate ternary complex through the catalytic rate constant k cat .
dNTP pool imbalance. The perturbation of components in the DNA precursor pool (denoted as Nuc p ) is influenced by aggregated effects of 5-FU stimulated TS inhibition, cellular nucleotide salvage pathway, and nucleotide feedback. Nuc p is represented by the variation of the total amount of deoxynucleotides at a given time point versus the corresponding control level, following the same methodology as [58]. The mathematical expression for Nuc p is determined as Eq (16).
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi In which x 0 and x denote the size of each constituent of dNTP pool before and after 5-FU treatment, respectively. The initial value of Nuc p is zero, implicating that prior to treatment there is no tendency towards unusual nucleotide biosynthesis. As the level of all dNTPs decrease or increase drastically from the baseline x 0 , the perturbation measure would increase accordingly.
The data used for model fitting is also converted to the measure of perturbation according to Eq (16). The time profiles of TS f and TS total , determined in the cellular model, are connected with changes of Nuc p by introducing a new variable called perturbation in TS inhibition (TS p ) into this submodel. The degree of TS inhibition is evaluated by the ratio of the number of free FdUMP binding sites (TS f , pmol�mg −1 ) and the total available TS binding sites (TS total , pmol�mg −1 ) [56]. TS p is quantified as the deviation of TS inhibition measured post-treatment (TS f /TS total ) from the normal portion of free TS levels (TS f0 /TS total0 ) (Eq (17)).
With initial condition of 1, the value of TS p varies within the range 0 to 1. TS p vanishes when TS nucleotide binding sites is fully occupied by FdUMP. We applied and modified the scheme previously developed by Setzer et al [58] to simulate the dNTP pool imbalance induced by 5-FU. The degree of dNTP pool imbalance is quantified based on the two assumptions: 1). The 5-FU cytotoxicity is directly related to dTTP and the alteration of dNTP is the direct and immediate product of TS inhibition; 2). Any disruption of dNTP supplies as a result of TS inhibition would impede the normal process of DNA replication. The kinetics of dNTP pool imbalance is decided by the two differential equations, one for variable Nuc p as discussed above and one for variable U. The cellular nucleotide salvage pathway is activated in the face of intra-cellular thymidine starvation caused by TS inhibition. Thus U is invented to describe the propensity of the cellular nucleotide salvage pathway to diminish the disturbance in the levels of constituents in the pool. The dynamics of U follow the pattern of a negative feedback loop (Eq (19)). As the first terms in Eqs (18) and (19) show, the rates of production of Nuc p and U both increase as the perturbation in TS inhibition are elevated, and finally converge to corresponding upper bounds(k 1 and k 10 ). Therefore, as no FdUMP is attached to TS, all the binding sites on TS are available for dUMP, and the fluctuation in dNTP would return to baseline level. In addition, the collective effects of TS p , Nuc p and U impose negative components to the rate of change of Nuc p and are represented by the third term in Eq (18) which is smaller for the smaller level of TS p and larger for the greater level of Nuc p and U. Besides, the second term in Eq (18), solely dependent on the level of Nuc p , describes the amelioration of disturbance by biochemical interactions among components in the dNTP pool. The kinetics of Nuc p and U are described by the differential equations(Eqs (18) and (19)).
In which γ dNTP is a hill parameter. All the other parameters are the unitless quantities. DSB induction. dNTP pool imbalance and cellular repair pathways are implicated in the production and elimination of DSBs, denoted as N DSB (t) in the model. The intactness of DNA precursors in the dNTP pool partially guarantees the success of DNA biosynthesis and DNA damage repair [32]. The perturbed dNTP pool caused by TS inhibition would result in double strand DNA break during DNA replication in the S phase [32,37]. In response to the aberrant DNA duplex structure, the cell cycle checkpoints are triggered to allow the repair machinery to restore the genomic stability; if the DNA repair fails, apoptosis pathway would be activated to clear the damaged cells [34].
The literature data shows that the induction of DSBs caused by 5-FU is delayed following the occurrence of dNTP pool imbalance. Such an assumption is also biologically reasonable since the two cellular events are occurred at different molecular levels and are related to cell cycle progression. In the absence of detrimental effects of dNTP pool imbalance, the generation of severe DNA lesions depends on the zero-order rate constant k 0 . The outflow term V max;HR �N DSB ðtÞ K m;HR þN DSB ðtÞ stands for the efficiency of normal DNA repair responses. Mathematically, the dynamics of DSB generation (N DSB (t)) for 0 � t � T d is governed by Eq (20) By contrast, the rate of formation of 5-FU-induced DSB depends on the past states of perturbation in dNTP pool (Nuc p ) and is modeled as a full sigmoid formula so that the amount of DNA fragmentation caused by 5-FU may converge to an upper limit V max,dNTP . The elimination of DSB is mediated not only by DNA repair response, but also by the availability of DNA precursors [59]. It has been suggested that the polymerization of DNA in the course of homologous recombination repair (HR) is obstructed by imbalanced nucleotide pools resulting from FdUMP binding to TS [60]. Therefore, the kinetics of HR efficiency decreases as the perturbation in the dNTP pool is elevated. For t > T d , the dynamics of DSB induction is governed by Eq 21.
In which the delaying effect of dNTP pool imbalance on N DSB (t) is modeled as the explicit delay parameter T d,DSB . γ DSB is used to adjust the shape of the curve.
5-FU-mediated disturbance in genomic stability (N DSB,deviation (t)) is qualified on the basis of kinetics of N DSB (t) and is incorporated into the TGI model. N DSB,deviation (t) measure deviation of the DSBs exclusively caused by dNTP pool contamination from the time-dependent baseline level, given by Eq (22) In which N DSB,deviation (t) is a unitless quantity. The time course of N DSB,0 (t) is simulated based on the placebo treatment [61]. The value of N DSB,deviation (t) is zero when no drug is administered and N DSB 0 ðtÞ is the number of DSB in the absence of any medical intervention and is varied with time as it is probable for other endogenous factors to induce DSBs to a lesser extent.

Tumor growth inhibition model
The TGI model is built on the basis of the indirect response model and life span model. It is assumed that the overall tumor cell population, denoted as T, comprises proliferative and damaged cells. The pharmaceutical response examined in this module is the change of tumor volume(cm 3 ). Accordingly, the tumor growth inhibition (TGI) model is composed of the differential equations governing the kinetics of proliferating cell population (P(t)) and damaged quiescent cell population (D(t)), as shown in Eqs (23) and (24).
Where T d,Tv is the mean life span of damaged cells. Cells D would be processing the damage for T d,Tv until they are eliminated from the entire tumor cell population. A 5 (t) is amount 5-FU anabolites at time t, computed in cellular model. IC 50 is the concentration of 5-FU anabolites that causes the 50% decrease in tumor proliferation rate. E max,damage is the maximal fractional factor of DSB-induced stimulatory effect on the production of damaged quiescent cells. EC 50,damage is the level of N DSB,deviation necessary to produce the half-maximal stimulatory effect. λ g is the tumor growth rate. λ d is the natural death rate. P max is the carrying capacity of tumor growth. γ Tv is a curve-adjusting parameter, whose value is fixed as 0.2. The initial value of P(t) is the last measurement of tumor volumes before administration starts, while the initial value of D(t) is zero.
The TGI model incorporates two mechanistic components, 5-FU anabolites and 5-FU induced DSBs to account for the inhibition of tumor proliferation by 5-FU treatment. Their individual effects on the pharmacological responses are mathematically formulated by full sigmoid E max,damage model, giving rise to E a (t) (Eq (27)) and E DSB (t) (Eq (28)) whose value is subject to the upper limit of 1. The smaller the value of IC 50 , the more responsive tumor cells would be to the suppressive effects imposed by 5-FU anabolites. Likewise, N DSB,deviation (t) is presumed to enhance the loss of proliferative cells with the rate given by E DSB (t) The maximum allowable value of E DSB (t) is E max,damage . E max,damage , EC 50,damage and IC 50 are drug-sensitivity parameters and modulate the potency of the chemotherapy. Two different scenarios associated with the TGI model are discussed in detail.

Tumor dynamics in the absence of drug intervention
The dynamics of proliferative cells in the absence of drug intervention is mediated by Eq (29).  ) is zero, meaning that no damaged cells emerge from normally cycling cell population. Taken together, the three parameters λ g , P max , λ d outline the tumor growth pattern in the absence of medical intervention.

Treatment responses
With 5-FU treatment, the tumor progression is assumed to be subject to 5-FU anabolites and DSBs induced by the perturbed dNTP pool. The effects of two quantities are modeled as hill functions and serve as system biomeasures driving tumor dynamics, as shown in Eq (24). The biochemical reactions mediated by 5-FU anabolites result in abnormality in DNA replication and repair, leading to the incomplete preparation for cell division. Such feature is reflected by the term E a (t) that is in inverse proportion to the concentration of 5-FU anabolites. Besides, as a contributory factor in 5-FU cytotoxicity, 5-FU-induced DSBs are responsible for the production of damaged cells with a rate E DSB (t). Therefore, in our model design, the increase of quiescent tumor cells is driven by the irreversible transition from proliferative tumor cells instead of their proliferation. The tumor cells that incur DSB lesions are not deemed to be eliminated from the entire tumor cell population instantaneously after the onset of DSB, and yet they would undergo progressive stages of damage before they die, which is biologically plausible since the processing of DSB is related with time-dependent signal cascades regulating cell cycle progression. To that end, the life span model, as reviewed in literature [62,63], is implemented to capture the elongated survival time of non-growing cells. As components of the life span model scheme, the delay differential equation (DDE) is incorporated in the TGI model and describes the delayed onset of the observed drug effects. T d,Tv denotes the mean transit time through all the damage stages for P(t). Therefore, in the case of 5-FU treatment, the alterations of variable A 5 (t), N DSB,deviation (t) directly change the dynamics of tumor growth. The damaged cells would emerge once 5-FU induced DSB is above the baseline value N DSB 0 ðtÞ. Between the two consecutive doses, N DSB, deviation (t) and A 5 (t) are expected to return to zero and the change in tumor volume would then become a function of λ g , P max , λ d , claiming the pattern of unperturbed tumor growth. Besides, the rate constant of tumor natural death and the baseline rate of irreversible conversion from the proliferative to quiescent cells are assumed to be equal toλ d .  (30):

Model implementation and Parameter estimation
where d ij ; i ¼ 1 . . . 4; j ¼ 1 . . . n i ; d ij 2 R L ij represents the time series data of the j th model variable(denoted as R ij ) whose kinetics are governed by the parameters in Θ i . The difference between each measurement and model prediction of R ij is described by the measurement noise ε ij . We assumed that ε ij are independent random variables and normally distributed with zero mean and the variance s 2 ij . Instead of keeping ε ij constant, we treated them as random variables whose marginal distributions can be learned from the data. f ij ðΘ i Þ 2 R L ij represents the simulated time course of j th variable using the parameter set Θ i . The choice of the prior distribution for each model parameter is listed in the Supporting Information S3 Text. It is assumed that all the parameters are random variables and are independent of each other. Thus, the joint posterior distribution of the model parameters in the parameter set Θ i , i = 1. . .4 can be determined by Eq (31).
where m i represents the number of parameters in Θ i . We used the robust adaptive Metropolis algorithm(RAM) [64] as the MCMC technique. The mean values of the estimated parameters together with the credible regions for each parameter are listed in Table 1. The histograms of the estimated posterior distributions for the model parameters are provided in the Supporting Information S3 Text.

Results
Model construction PK model. The PK model is composed of central and peripheral compartments. The parameters in the PK model are determined by fitting the model to the plasma concentration data measured on animal models injected with a single 100 mg/kg 5-FU dose [65]. To maintain unit consistency throughout the integrated model, the unit of 5-FU plasma concentration is converted to a molarity concentration unit pmol/mL from μg/mL via the molar mass of 5-FU (130.077 g/mol). The time course of 5-FU in the peripheral compartment is estimated from the PK data. The parameter estimates are shown in Table 1 and simulated concentration-time curves for compartments 1 and 2 are shown in Fig 2. From the result of parameter estimation, the estimated volume of the central compartment is smaller than that of the peripheral compartment, which agrees with the estimation results shown in the previous PK/PD study [20].

Computational analysis of 5-fluorouracil anti-tumor activity
The elimination phase of the peripheral compartment is slower than that of the central compartment. Cellular model. The cellular model is composed of 5-FU concentration in interstitial fluid (Eq (4)), intra-cellular 5-FU(Eq (5)) formation of 5-FU anabolites (Eq (6), insertion of FUTP, FdUTP into RNA and DNA (Eqs (7) and (8)), binding of FdUMP to TS (Eq (9)) and dUMP (Eq (14)), perturbed dNTP pool caused by reduced TS catalytic activity (Eq (18)) and DSB induction (Eq (19)), whose trajectories are expressed by the solutions of the differential equations. The parameters in the cellular model are estimated by fitting to data measured on colon xenograft tumors. The literature data sources are enumerated in Table 2 and the parameter estimates are provided in Table 1. The initial values of the variables are derived from literature data, that is the content of species at 0 h upon exposure to 5-FU. Instead of fitting to literature deta, the time course of A 3 (t) is estimated dependent on the fluxes between compartment 3 and the neighboring measurable quantities.
As shown in the model simulations (Fig 3C), the accumulation of 5-FU anabolites begins within minutes upon exposure to 5-FU, subsequent to the transmission of 5-FU from interstitial fluid to intra-cellular space. The model is able to capture the sharp decrease of non FdUMP bound TS enzyme in parallel to a rapid increase in dUMP, indicating pronounced inhibition of TS activity. In Fig 3F and 3G, it can be observed that the peak of dUMP amount and the lowest level of free TS enzyme coincide at roughly 5 hours after 5-FU administration, after which the decline in dUMP is accompanied by the restoration of free TS. In terms of 5-FU incorporation into the DNA, as demonstrated in both model simulation (Fig 3D and 3E) and literature data, the fraction of FU-RNA is higher than that of FU-DNA. Such a phenomenon may be related to the fact that RNA replication happens all the time in the cell while DNA replication only happens in the S phase. In addition, the duration of incorporation of 5-FU in   RNA and DNA is longer than that of other substances, which can be explained by the involvement of cell cycle kinetics in regulating the genomic activity. Fig 3I depicting the formation of DNA lesions caused by 5-FU shows that in the early stages following the 5-FU administration, DSBs are maintained at a certain level before their count noticeably increase after approximately 40 hours and ultimately reaches an equilibrium state. In comparison with the timing of onset of other species, the appearance of 5-FU-induced DNA fragmentation is delayed, because 5FU can only exert anti-tumor effects when DNA replication takes place. By visually comparing Fig 3G and 3H, we also noticed that the maximum values of dNTP pool imbalance and 5-FU-induced DSB are reached asynchronously, in terms with the observations shown in the experimental studies focusing on the impacts of dNTP pools on the generation of lethal DNA damage and cell fate [37]. However, the time elapsed between the peak of dNTP pool imbalance and the generation of 5-FU-initiated DSB may be different if the model is fit to distinct colon cell lines with differential growth pattern and initial cell cycle phase. The histograms of the cellular model parameters are given in the Supporting Information S3 Text. We classified the parameters into identifiable and unidentifiable groups according to the compactness of the supports and shapes of their distributions. As shown in Fig B and Table 3 to probe the alteration in the dose-related parameters as the changes of doses and injection schedules. The parameters in PK and cellular model are kept the same as the ones obtained from colon xenograft tumors given a single 100 mg/kg 5-FU. The initial values of 5-FU in central plasma are varied according to the dosing schemes. The PK model is reloaded at the frequency of dose administration. The parameters describing the kinetics of control groups are estimated for all the dose regimes whose data are from various sources in the literature, keeping the same value for γ Tv . The Fig 4 demonstrates that the model simulations using the mean parameter estimates shown in Table 4 are consistent with the tumor volume dynamics in literature data. Moreover, due to the lack of repeated experimental measurements from the plasma-treated (PBS) control groups, the zero net growth rate around the maximum level is assumed, indicating that the model parameterization could not explain tumor growth above the equilibrium. As the drug effect wears off, the tumor volume would eventually approach an asymptote equal to P max , the value of which is directly estimated from the time series data with the limited duration. For dosage regimens (1), (2), (3), (5), the parameters in the TGI model are estimated by fitting the literature data, forming models (1), (2), (3), (5). Following this, the parameter estimates for the TGI model 3 are reused on protocols (4), (6), and (7) except for the drug-irrelevant parameters, considering that they have the (2) 50 mg/kg, daily for 4 days Colon 26 [72] (3) 20 mg/kg, twice a week for 3 weeks HT 29 [73] (4) 20 mg/kg, three times a week for 3 weeks HT 29 [74] (5) 5 mg/kg, daily for 9 days SW620 [75] (6) 20 mg/kg, twice a week for 3 weeks SW620 [76] (7) 20 mg/kg, every other day for 2 weeks SW620 [77] (3) is also applied to models (8) and (9). On the other hand, the estimation of TGI parameters except for T d,Tv for the protocols (8) and (9) are carried out separately due to their unique dosing structures and type of cell line. Besides, the parameters governing normal tumor growth pattern (i.e. λ d , P max , λ g ) are estimated for each protocol since the experimental data for each protocol is from a different study in the literature, and is subject to the characteristics of the cell lines and progression phases of the xenograft tumors. The control group data is used to characterize the parameters λ d , P max , λ g . The literature tumor volume data are normalized to the respective starting values. It is also assumed that for each model, the first doses are injected at 0 h after the treatment course begins. The values of parameters shown in Table 4 also agree with the extent of adverse effects of the studied dosage regimens. Among all the nine models, the estimated value of IC 50 for model (1) is the largest. It is partially due to more anabolites that are synthesized from a higher level of drug exposure to tumor cells. The DSB-related parameters E max , EC 50,damage of model (2) are the highest, which is caused by the acute accumulation of DNA lesions as a result of daily administration. As for the mild dose schedule featuring models (1) and (3), the amount of single dose plays a more dominant role in determining the levels of E max and EC 50,damage , which is indicated by the larger value of the two model parameters for injection of 100 mg/kg 5-FU than that of 20 mg/kg 5-FU. The values of those growth-related parameters are subject to the characteristics of the cell lines and progression phases of the specific xenograft tumors. Overall, the parameter estimates of 9 xenograft tumor models examined in this study fall in the same order of magnitude. Besides, the computational analysis of 5-FU anabolites-effect relationship (E a (t)) and DSB-effect relationship (E DSB (t)) are documented in Supporting Information S1 Text. Fig D-Fig L in Supporting Information S3 Text show the histograms for the estimated TGI model parameters. The identifiable and unidentifiable parameters for the TGI models (1)  Table 4. Red line represents the time course of the control group; blue line represents the time course of colon tumor growth treated with corresponding dosage regimens shown in Table 3; triangle, literature data for the control group; circle, literature data for the treated group; error bar, standard errors; asterisk on x-axis represents injection time. The shaded areas in each panel represent the 0.7 credible intervals.  -(9) are listed on the figure captions in Supporting Information S3 Text. We conducted a global sensitivity analysis where we first used the Morris method to screen the parameters that have more influence on the model outputs. Then the variance-based global sensitivity analysis, often referred to as Sobol method, is conducted based on the selected parameters. The Morris method revelaed that the changes in the parameters K m , K m,efflux , K m,54 , K m,influx , P max , V max, 75 , V max,influx , V max , λ d , λ g , k 03 , k 56 , k 59 , k 95 , k 65 created considerable perturbations on the model output and these were chosen to be analyzed using the Sobol analysis. The detailed discussion of the analysis and the results are provided in Supporting Information S2 Text. We further used the Sobol global sensitivity to identify the degree of the changes in the model outputs produced by the joint interactions of the input parameters. The first order and total order Sobol index of each parameter with the tumor volume at 5th, 25th and 45 day after treatment begins are given in Supporting Information S2 Text. The results show that model outputs at 45 th day are more sensitive to λ g , P max , V max,efflux , and V max,influx than the other input factors, which is consistent with the result of preliminary screening by Morris method.

Quantitative evaluation of treatment efficacy using established TGI models
We further examined extensively the effect of dosage regimens on cytotoxic effects of 5-FU by looking at the area under the curves of E a (t) and E DSB (t) (denoted as AUC e,a and AUC e,damage respectively), the area under the time profiles of A 5 (t) and N DSB,deviation (denoted as AUC a and AUC DSB respectively), the trajectories of E a (t) and E DSB (t) (shown in S1 Text), and tumor responses (Fig 4). The numerical results of AUC e,a , AUC e,damage , AUC DSB , and AUC a are shown in Table 5. AUC e,a and AUC e,damage are used as the indicators of the cumulative cytotoxic effect of the drug and can be used to predict cumulative tumor responses. Following the model design, AUC e,damage increases with the accumulation of drug exposure, while AUC e,a experiences the opposite trend. For the protocols (1), (3), (4), (6), (7), their accumulative doses over a week are similar, so is the accumulative effect of E a (t) and E DSB (t). But by close comparison of the individual AUC e,damage , the DNA damage caused by 100 mg/kg weekly injection exerts more significant anti-tumor effects among those regimens. And the treatment with 20 mg/kg 5-FU every other day (model (7)) stands out with the smallest AUC e,a . Furthermore, through the comparison of the AUC measures calculated for the protocol (2) and protocol (8), we can observe that the moderate dose coupled with the intense dose frequency can lead to striking anti-tumor effects. The AUC DSB and AUC a corresponding to daily 50 mg/kg 5-FU treatment are the extremest in comparison with those of other schemes. Likewise, the associated trajectory of the relative tumor volume change also demonstrates that daily administration of 50 mg/ kg 5-FU for 4 times shows more noticeable cytotoxic effects and causes the tumor shrinkage  within the first week after the treatment begins. Such phenomenon can be explained by the observation that the intracellular drug action mediators with extensive amounts resulted from frequent injections are present and continue to take effect over a protracted period of time.
The simulation results of models (3), (4), (7) and the corresponding calculated AUC e show that given the same size of dose and period of treatment, the more frequent drug injection results in the more eminent drug efficacy. It is also worth looking at the time profiles of E DSB (t) and E a (t) for model (1) to assess the combinations of high dose and less frequent schedule. It can be observed that while the drug effects would fade away during the dose intervals, high doses with a weekly schedule can still yield drastic tumor suppression. Taken together, 5-FU suppresses cell proliferation of colon cancer in a time-and dosedependent manner and the accumulative PK exposure can improve the overall therapeutic efficacy. The treatment with daily injection of 50 mg/kg 5-FU four times exhibits the most detrimental effect on tumor burden. A high-dose 5-FU taken over long period may produce significant anti-tumor effects. Regarding the same doses given over the same treatment duration at specific intervals, the more intense schedule may increase the likelihood of exposing vulnerable tumor cells to 5-FU and yields a more significant anti-tumor effect.

Simulation analysis of tumor resistance towards 5-FU exposure
In this section, we investigate various scenarios on the mechanisms underlying the resistance of tumor responses towards 5-FU treatment using numerical simulations. We evaluated simultaneous changes in physiologically-related parameters that would generate a resistant behavior. We summarized the parameter variations we used to identify possible resistance mechanisms on Table 6, which also lists how these variable changes reflect the biologically relevant resistance mechanisms studied in the literature. The last three entries on the table show our novel hypothesis on other possible causes that might lead to treatment resistance. Considering that resistance may be caused by interwoven cellular pathways, we simulated the simultaneous changes of certain factors to examine the aggregate effects. The Monte Carlo method is applied to generate N(= 500) parameter sets and time courses of tumor volume change are simulated for each parameter set. Each parameter set comprises the concerned parameters with newly sampled values from their associated sampling intervals(Eq (32)) and the Table 6. The range of variation for each parameter and their corresponding biological implication.  [56]; metabolism may predict tumor sensitivity to 5-FU therapy [40] K m, 54 2 fold increase -Decreased accumulation of anabolites G 0 10 fold increase 10 fold decrease dUMP accumulation [56] k cat 10 fold increase 10 fold decrease dUMP accumulation k 08 10 fold increase 10 fold decrease dUMP accumulation k 95,complex -10 fold decrease decreased enzyme affinity for FdUMP [56,79] k 59,complex 10 fold increase -Decreased stability of ternary complex [56,79] k 09,complex 10 fold increase -Decreased stability of ternary complex [56,79] k 5,dNTP 10 fold decrease The effects of salvage pathway on lessoning dNTP pool imbalance [80].
T d,DSB 2 fold increase -Delayed 5-FU-induced DNA fragmentation and cytotoxicity in the resistant cell lines [79].  (8) and (9). The resistant tumor response is assumed to be the simulation of the TGI model (9). The parameters we experimented on in this section are the ones related with 5-FU drug actions and resistance mechanisms, which are V max,54 , K m,54 , G 0 , k cat , k 08 , k 59 Where p l i and p u i are the lower bound and upper bound of parameter p i . rand represents a random number drawn from the standard uniform distribution. The ranges of variation for each parameter are shown in Table 6.
The departure of each tumor growth time course simulated using a new set of parameters from the time course for resistant LS174T tumor cells (model (9)) is calculated as the sum of the Euclidean distances between the coupled points on the two time course trajectories. The time course closest to the insensitive / resistant tumor response is the one with the smallest departure. The mean time course for a certain case is calculated by averaging all the computed time courses. Maximal and minimal time courses among all the simulation results are determined by calculating the ℓ 2 -norm of Euclidean distances between the sample points along certain path and x-axis. Table 7 demonstrates parameter variations that can give rise to tumor response closest to the resistant LS174T xenograft tumor(model (9)) for all the seven cases A-G.
The results that are given in

PLOS COMPUTATIONAL BIOLOGY
Computational analysis of 5-fluorouracil anti-tumor activity significant effect by plummeting N DSB to its baseline. The cases (D) and (G) demonstrate similar resistance responses in terms of the tumor volume and magnitude of 5-FU induced DSBs.

Generalizability of the mechanism-based model
The generalizability and accuracy of our proposed mechanism-based model on the prediction of tumor responses given a specific dose and dosage regimen are evaluated in this section. We applied different dosing schedules to the TGI model developed above to see how well it can capture the literature data for the target dosing schedules. The values of efficacy-related parameters (i.e. E max , EC 50,damage , IC 50 , T d,Tv ) are kept the same for models (3), (4), (6), (7), because these models have the same amount of single dose, implying that the levels of drug-related intermediates are expected to be the same for the single administration. Because the cell line differences have been addressed through control groups, as shown in Table 4, the distinct dose schedules should be the sole drug-regulated factor that drives the differences in the outputs of these models. Therefore, the efficacy-related parameters (i.e. E max , EC 50,damage , IC 50 , T d,Tv ,) of models (4), (6) and (7) are kept equivalent to those of model (3). In order to facilitate the assessment of model generalizability, we organized the TGI models except for models (1), (2), (8) and (9) into three different groups according to the rules defined as follows: (A) models with the same type of xenograft tumor and the same single dose but with different frequencies; (B) models with a different type of xenograft tumor under the same protocols; (C) models with a different type of xenograft tumor under the same amount of dose but different frequencies. Models (3) and (4) fall under scenario A. Models (3) and (6) satisfy the condition of scenario B. The last scenario C involves models (3), (4), (6), (7). Fig 4D indicates that the model 4 simulations for protocol (4) given in Table 3 moderately captures the data, which suggests that the model is structurally stable in the face of changes in the dose frequency. Regarding scenario B, the selections of protocols (3) and (6) from Table 3 show that the differences in treatment responses are caused by the distinction in the types of xenograft tumor. Fig 4C and 4F show that the model simulations agree well with literature data, indicating that our model can adapt to the different tumor types merely through adjusting the parameters for normal tumor growth. Lastly, the experimental measurements of model (3), (4), (6), and (7) suggest that changing dose schedules and types of cell lines influence tumor responses to multiple doses, which parallels the simulations of these four models. The model outcomes of (3), (6) and (7) successfully capture the experimental dynamics, whereas model (4) has a moderate fit to the experimental data. From the generalizability standpoint, model (4) that shares the same efficacy-related parameters was able to capture the general trend in the observed data. Based on the outcomes of evaluation discussed above, it can be concluded that the proposed model constructed under a certain dosing scheme is capable of predicting the tumor responses towards a similar dosing scheme or the one with the same amount of dose but different administration frequency, as long as the control tumor growth has been captured by the model with no drug input.
We utilized our models to computationally experiment with dose regimens that are different than the ones given in Table 3 to further examine the anticancer efficacy of 5-FU treatment. We developed four hypothetical regimens with different combinations of doses, dosing frequencies and treatment durations than the ones we used from the literature. The regimens used for these computational experiments are (A) 100 mg/kg on day 1 of weeks 1 to 6 of an eight-week cycle, for a total of 4 cycles; (B) 50 mg/kg twice a week of weeks 1 to 6 of an eightweek cycle, for a total of 4 cycles; (C) 100 mg/kg on day 1 of weeks 1 to 3 followed by 50 mg/kg twice a week of weeks 4 to 6, for a total of 4 cycles (32 weeks); (D) 50 mg/kg twice a week of weeks 1 to 3 followed by 100 mg/kg on day 1 of weeks 4 to 6, for a total of 4 cycles (32 weeks).
The simulation results are shown in Fig 6. Despite having the same total amount of doses over a week, dosage regimens (A) and (B) result in notably different treatment outcomes due to the difference in administration times. The simulation result in Fig 6B shows that within the first cycle, the consecutive twice-a-week injection of 50 mg/kg 5-FU causes a significant tumor shrinkage while the tumor still grows under weekly injection of 100 mg/kg as shown in Fig 6A.  Fig 6C and 6D demonstrate that the order of administration of different dosage regimens makes a difference for the treatment options comprising mixed doses. Injection of 50 mg/kg 5-FU twice a week for three weeks following weekly injection of 100 mg/kg 5-FU can inhibit the tumor growth more than the regimen with the opposite order. Computationally, all of the dose regimens presented here show more favorable treatment outcomes compared to the TGI model (1) with the dose regimen 100mg/kg administered weekly for four weeks. Comparison between TGI model(1) and dosing regimen (3) on relative tumor change between three weeks and seven weeks after treatment reveals that shortening the intervals between consecutive lower dosages might yield better outcomes of cytotoxicity. It should be noted that these dose regimens are developed to conduct computational experiments on treatment outcomes and it is necessary to experimentally identify the toxicity profiles and the maximum tolerable doses for the evaluation of dosing schedules before they can be used in a clinical setting. It can be concluded from the computational experiments that applying dosage regimens with high potency after the ones with less potency may result in clinically favorable treatment outcomes.

Discussion
Mathematical modeling is a powerful tool for capturing the target pharmacology and making predictions on tumor response by translating current understanding towards underlying mechanisms into rigorous mathematical terms. Our physiology-based model allows for the analysis of the dynamics of critical intermediate components related to 5-FU effective drug actions on different time scales and across physical domains. Our cellular model quantitatively reflected the drug effects and related the mechanisms of TS inhibition with disturbance in the dNTP pool and induction of DSBs. The tumor heterogeneity was integrated into the TGI model. The combination of inhibition of proliferating progenitor cells and production of damaged quiescent cells characterized the dynamics of suppressed tumor growth. We converted the time profiles of the 5-FU anabolites and 5-FU-induced DSBs into nonlinear measures of drug effect on inhibiting tumor growth rate and simulating the loss of proliferative cells to provide a realistic prediction of the kinetics of tumor progression under 5-FU treatment. The integrated computational model has the potential to promote the identification of contributions of the physiological intermediates on the clinically favorable responses and potentially enable the precise prediction of medical treatment by antagonistic agents.
Like the other mathematical models, our model is inevitably influenced by a certain level of uncertainty. The uncertainty may stem from the insufficient measurements, parameter estimation, and alternative model structures, which propagate and aggregate as the model progresses. Furthermore, the solution of the dynamic model could be distracted by a local minimum, in the sense that the model is inherently un-convex. This means that we cannot solely depend on the point estimates provided by the optimization scheme involving the minimization of the objective function. Instead, gaining the knowledge of the joint posterior distribution of model parameters is beneficial in terms of recognizing model identifiability. To that end, we implemented the sampling-based Bayesian inference method through the MCMC technique. The model identifiability was examined by plotting the histograms for each parameter in the sample set. In addition, although the time course measurements for system variables are collected from different literature, we attempt to constrain the scope of data to in vivo or in vitro measurements under 100 mg/kg 5-FU treatment.
Examination of the combined effect of doses, dosage regimens and duration of treatment via the established model helps us seek the concealed relationship among the treatment protocols. When simulating different dosage regimens, we adapt to the lack of time-series data of fluorinated compounds under other amounts of drug dose than 100 mg/kg by using the cellular model obtained using 100 mg/kg data collections. This may be due to the fact that the 5-FU anabolites possess a rapid rate of transformation, and capturing them under smaller doses to a precise degree may run into technical difficulties. Therefore, it is unrealistic to reevaluate the parameters of the cellular model for each dosing regimen. We used the cellular model determined by fitting to 100mg/kg literature data to simulate all the aforementioned dose and dose schedules. From the model simulations under multiple dosage regimens, we found that a single high dose outperforms the divided dose in terms of drug efficacy of causing tumor shrinkage. Besides, the more intense schedules are able to yield higher therapeutic efficacy compared to the protocols with the same doses and treatment course but less frequent schedules. Taken together, our model can provide an accurate prediction of pharmacological responses to different dosage regimens. Through such investigation, based on our comprehensive model, it has the potential to advance the quantification of chemotherapy potency and to facilitate the development of efficacy optimization.
We also extended the application of the model to the analysis of tumor resistance responses to the cellular defense mechanisms attributed to tumor insensitivity towards 5-FU, influencing the final treatment outcome. For instance, deoxyuridine triphosphate nucleotidohydrolase (dUTPase), an enzyme catalyzing the degradation of FdUTP to FdUMP, plays a contributory role in enhancing tumor resistance. The experimental study suggests that the cellular (F)dUTP accumulation resulting from dUTPase activity may be related to the insensitivity of tumor to fluoropyrimidine-induced cytotoxicity [79]. Moreover, based on the experiments conducted on mice bearing murine mammary carcinoma, the researchers have found that the incorporation into RNA can be blocked by deoxyuridine (dUrd), rescuing the tumor cells from the lethal toxicity of FUra [81]. Homologous recombination repair (HR) can protect the tumor cells against exogenous insults. DSBs can be repaired by HR with high fidelity [82]. The expression level of the HR proteins is also related to the tumor resistance to the 5-FU anti-tumor effects [60]. In addition, it has been proposed that a slow-down in cell cycle and the reduced expression of CDK2 protein can prevent the incorporation of 5-FU metabolites into DNA [56]. However, all the cellular mechanisms described are more or less relevant to 5-FU -induced TS inhibition since TS is a key player in cellular nucleotide metabolism. Therefore, we conducted the analysis of tumor resistance responses around 5-FU induced TS inhibition and its downstream events, such as the dNTP pool imbalance and DSB induction. As a result of our computational analysis, the decreased accumulation of 5FU anabolites, variations of dUMP from baseline, reduced production of TS-FdUMP complex, and eased pressure from 5-FU induced DSB may bring about the tumor resistance responses.
Taken together, our integrated model forms a nested multi-scale structure and quantitatively connects the effects of 5-FU on tumor growth with genome instability and PK-dependent physiological turnover processes, facilitating simulation studies and resistance analysis. The PK, cellular, and tumor growth inhibition models were arranged in accordance with the mechanistic sequence, providing a realistic description of physiological events that ultimately decide the tumor fate. By fitting to literature data measured from preclinical animal models under different dosage regimens, the proposed model showed robust predictive performance and can simulate the PK and PD profiles under different protocols, reflecting that the model has the potential of being applied in the protocol designing phase. Our proposed semi-mechanistic PK/PD approach provides insight towards translation from xenograft tumor study to clinical efficacy analysis in a computational manner. The model simulations can also provide a tool for researchers to compare and evaluate the impacts of various possible dose and dose frequencies on the improvement of drug efficacy, driving the selection of desirable chemotherapy treatment.
Supporting information S1 Text. Quantitative analysis of biomeasures E a (t) and E DSB (t) in TGI model. The file details the computational analysis of 5-FU anabolites-effect relationship(E a (t)) and DSB-effect relationship (E DSB (t)).